Comparison of the Transcriptomic Atlas of the Human Retina and Mouse Retina

Setup

Load libraries

library(ggplot2)
library(tidyr)
library(dplyr)
library(Matrix)
library(Seurat)
library(cowplot)
library(patchwork)

Process Human Data

import_remote_data <- function(file_url, type = "table", header = FALSE) {
  con <- gzcon(url(file_url))
  txt <- readLines(con)
  if (type == "MM") { return (readMM(textConnection(txt))) }
  if (type == "table") { return (read.table(textConnection(txt), header = header)) }
}

count_matrix_URL <- "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE137nnn/GSE137537/suppl/GSE137537_counts.mtx.gz"
gene_names_URL <- "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE137nnn/GSE137537/suppl/GSE137537_gene_names.txt.gz"
sample_annotations_URL <- "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE137nnn/GSE137537/suppl/GSE137537_sample_annotations.tsv.gz"

count_matrix <- as.matrix(import_remote_data(count_matrix_URL, type = "MM"))
gene_names <- import_remote_data(gene_names_URL, type = "table")
sample_annotations <- import_remote_data(sample_annotations_URL, type = "table", header = TRUE)

rownames(count_matrix) <- tolower(gene_names[,1])
colnames(count_matrix) <- tolower(sample_annotations[,1])

human_ret_seurat <- CreateSeuratObject(counts = count_matrix, 
                                       meta.data = sample_annotations, 
                                       project = "human_ret", 
                                       min.cells = 3, 
                                       min.features = 200)

Some cells in meta.data not present in provided counts matrixNon-unique features (rownames) present in the input matrix, making uniqueFeature names cannot have underscores ('_'), replacing with dashes ('-')

human_ret_seurat

An object of class Seurat 
19712 features across 20091 samples within 1 assay 
Active assay: RNA (19712 features, 0 variable features)

Process Mouse Data

mouse.data <- Read10X(data.dir = "filtered_feature_bc_matrix")
dimnames(mouse.data)[[1]] <- tolower(dimnames(mouse.data)[[1]])
dimnames(mouse.data)[[2]] <- tolower(dimnames(mouse.data)[[2]])
mouse_ret_seurat <- CreateSeuratObject(counts = mouse.data, 
                                       project = "mouse_ret", 
                                       min.cells = 3, 
                                       min.features = 200)
mouse_ret_seurat

An object of class Seurat 
16424 features across 4510 samples within 1 assay 
Active assay: RNA (16424 features, 0 variable features)

Combine

ret.list <- list(human = human_ret_seurat, mouse = mouse_ret_seurat)

ret.list <- lapply(X = ret.list, FUN = function(x) {
    x <- NormalizeData(x, verbose = FALSE)
    x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000, verbose = FALSE)
})

Integration

ret.anchors <- FindIntegrationAnchors(object.list = ret.list, dims = 1:50,  anchor.features = 1000)
ret.combined <- IntegrateData(anchorset = ret.anchors, dims = 1:50)

Integrated Analysis

DefaultAssay(ret.combined) <- "integrated"

# Run the standard workflow for visualization and clustering
ret.combined <- ScaleData(ret.combined, verbose = FALSE)
ret.combined <- RunPCA(ret.combined, npcs = 50, verbose = FALSE)
# t-SNE and Clustering
ret.combined <- RunUMAP(ret.combined, reduction = "pca", dims = 1:35)
ret.combined <- RunTSNE(ret.combined, reduction = "pca", dims = 1:35)
ret.combined <- FindNeighbors(ret.combined, reduction = "pca", dims = 1:35)
ret.combined <- FindClusters(ret.combined, resolution = 0.05)

UMAP Visualization

# Visualization
p1 <- DimPlot(ret.combined, reduction = "umap", group.by = "orig.ident")
p2 <- DimPlot(ret.combined, reduction = "umap", label = TRUE)
plot_grid(p1, p2)

DimPlot(ret.combined, reduction = "umap", split.by = "orig.ident")

TSNE Visualization

p1 <- DimPlot(ret.combined, reduction = "tsne", group.by = "orig.ident")
p2 <- DimPlot(ret.combined, reduction = "tsne", label = TRUE)
plot_grid(p1, p2)

DimPlot(ret.combined, reduction = "tsne", split.by = "orig.ident")

Identify Conserved Markers

max_cluster <- 8
marker.list <- as.list(c(0:max_cluster))
marker.list <- lapply(X = marker.list, FUN = function(x) {
  x <- FindConservedMarkers(ret.combined, ident.1 = x, grouping.var = "orig.ident", verbose = TRUE)
})

DefaultAssay(ret.combined) <- "RNA"
plots <- as.list(1:(max_cluster+1))
plots <- lapply(X = plots, FUN = function(x) {
  x <- FeaturePlot(ret.combined, features = rownames(marker.list[[x]])[1:6], min.cutoff = "q9")
})

for (i in seq_along(plots)) {
  t <- sprintf("Cluster %d", (i-1))
  p <- plots[[i]] + plot_annotation(title = t)
  print(p)
  #print(sprintf("Cluster %d", (i-1))
  #plot(plots[[i]])
}